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ABSTRACT 

Aims. We present a general method to calculate radiative transfer including scattering in the continuum as 
well as in lines in spherically symmetric systems that are influenced by the effects of general relativity (GR). 
We utilize a comoving wavelength ansatz that allows to resolve spectral lines throughout the atmosphere. 
Methods. The used numerical solution is an operator splitting (OS) technique that uses a characteristic 
formal solution. The bending of photon paths and the wavelength shifts due to the effects of GR are fully 
taken into account, as is the treatment of image generation in a curved spacetime. 

Results. We describe the algorithm we use and demonstrate the effects of GR on the radiative transport 
of a two level atom line in a neutron star like atmosphere for various combinations of continuous and line 
scattering coefficients. In addition, we present grey continuum models and discuss the effects of different 
scattering albedos on the emergent spectra and the determination of effective temperatures and radii of 
neutron star atmospheres. 
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1. Introduction 

In the last few decades much effort has been put into 
the modeling of radiative transfer in relativistically 
moving atmospheres such as novae and supernovae. 
One of the state of the art techniques to solve this 
problem is the operator splitting (OS) method. This 
method has been successfully used to solve radia- 
tive transfer problems with scattering and complete 
treatment of non-local thermodynamic equilibrium 
(NLTE) effects. 

However, these sophisticated methods for treat- 
ing radiative transfer have not been used in a gen- 
eral relativistic environment, although the form of 
the GR equation of transfer does, basically, not dif- 
fer from the special relativistic version and the OS 
method can be applied to such problems. 

There have been several successful attempts to 
model the emergent spectra of general relativistic 
systems such as neutron stars. These models also ac- 
count for magnetic fields or different surface temper- 



atures. However, the radiative transfer in these cases 
gen erally solves classical plane parallel problems. 
See lZavlin & Pavlovl d2002l) for a review. 

Here we present an OS method of solv- 
ing GR radiative transfer problems in spheri- 
cal (Schwarzschild) geometry. Other authors have 
solved the radiative transfer problem in GR. For in- 
stance Schinder & Bludma nf (I1989I) solve the mo- 
ment equations of radiative transfer with a variable 
Eddington factor method. The advantage of OS is not 
only that it does not depend on closure conditions but 
can also solve magneto-optical transfer, what could 
prove to be important as far as neutron stars are con- 
cerned; 

Zan e et all (Qj96) developed a characteristics 
method to solve general relativistic radiative trans- 
port problems. They utilize the constants of mo- 
tion for the description of photon orbits that arise 
due to the Killing vectors of the spherically sym- 
metric spacetime and use the analytical connection 
of an affine parameter with the radial coordinate as 
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well as choosing the momentum variables along the 
characteristics to be constant to formulate the radia- 
tive transport equation. Although this simplifies the 
equation in the case of flows such as accretion onto 
black holes and neutron stars, the lack of a comov- 
ing wavelength description forces the use of a large 
number of characteristics to resolve the (in this case) 
angular dependent absorption coefficients and a huge 
number of wavelength points to resolve the shift of 
spectral lines through the atmosphere. The latter is 
necessary even in the case of static general relativis- 
tic atmospheres, e.g., in neutron stars. 

We avoid these problems by using a comoving 
wavelength coordinate which explicitly accounts for 
the coupling of different wavelengths. This ansatz is 
the important core part of this work. Since in order to 
calculate detailed spectra you need to resolve spec- 
tral lines throughout the atmosphere and you cannot 
afford to use a wavelength description that depends 
on the layer of the atmosphere that you are at. Since 
in order to perform NLTE calculations you would 
have to add a significant number of wavelengths to 
your computation to resolve the spectral line at hand 
in every layer with the desired quality. 

In the following we present calculations of gen- 
eral relativistic radiative line and continuum trans- 
fer with a complete treatment of scattering. In order 
to demonstrate the functionality of the method, we 
chose a simple test case similar to a neutron star at- 
mosphere. 

2. Radiative Transfer 

Lindquist found the equatio n for radiative t rans- 
fer for a comoving metric dLindquistl 119661) . He 
used the photon distribution function as variable 
to describe the radiation fiel d. The wo r ks of 
ISchinder & Bludmanl d!989l) and IZane et all d 19961) 
follow this ansatz. 

However, we want to utilize a description of the 
radiation via the specific intensity that is suitable for 
our method of solution. 

Rather than using the general metric used by 
Lindquist, we neglect the effects of the atmosphere 
on the metric and use the Schwarzschild solution: 
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Since the atmospheres that are influenced by GR ef- 
fects typically have a small mass compared to the 
parent object, this simplification is well justified. 
Furthermore, it is possible to calculate the metric 
coefficients for spherical symmetry by integrating 
the Tolman-Oppenheimer- Volkov equations and thus 
avoiding this approximation if desired. 

The equation of radiative transfer in the 
Schwarzschild metric is then found in its character- 
istic form: 
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This development is equivalent to the work of 
Lindquist and delivers no new physical insight be- 
sides being in the same form as the spherically 
symmetric spec ial relativistic equation of transfer 
dMihalasl Il980l) . Hence modern operator splitting 
techniques to solve radiative transfer are applicable 
to the problem. 

However, there is a fundamental difference from 
the special relativistic equation of transfer. The co- 
efficient a A does not change sign in monotonic flows 
that describe, e.g., supernova and nova atmospheres. 
In the GR case, this coefficient is linear in jx and 
hence the sign of a A will be different for ingoing 
(yU < 0) and outgoing (p > 0) photons. a A cou- 
ples different wavelengths and determines how the 
spectral features shift due to the flow or the gravi- 
tational field. In the case of a supernova atmosphere 
the different parts of the atmosphere all move away 
from each other so that the direction of the wave- 
length shift along a ray is always the same, hence 
a A has the same sign along the ray. In a gravitational 
field, ingoing photons will be blueshifted and outgo- 
ing photons will be redshifted and hence, the sign of 
a A changes. This presents a difficulty as the direc- 
tion of flow of information in the wavelength space 
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is reversed along a ray and the transfer equation is no 
longer an initial value problem in wavelength space, 
but a 2-point boundary value problem. 

Mihalas already realized this in Mihalas ( 1980) 
and outlined a simple ray by ray formal solution 
to this problem. However, this solution is of lit- 
tle use for the construction of an approximate A- 
operator for an ALI-iteration. To solve this prob- 
lem we use the OS meth od described in detail in 
iBaron & Hauschildj (120041) . 

The treatment of arbitrarily flows of informa- 
tion in the wavelengthspace means that every spatial 
point at all wavelengths can influence the intensity 
at a given spatial point at a given wavelength. This 
forces you to use a matrix notation for the formal so- 
lution where you have to implement proper bound- 
ary conditions gouverned by ax at every spatial point 
for all wavelengths to ensure a locally stable upwind 
scheme for the wavelength discretization. 

The GR case is a rather simple application to this 
method since the flow of information in the wave- 
lengthspace just changes once on a ray, namely at the 
point of tangency or in the case of core intersecting 
rays at the innermost layer. 

To make this method work with the GR case 
you have to calculate ax for every wavelength on all 
points of a given characteristic and determine the ap- 
propriate discretization of the derivative at this point. 
Furthermore you have to know the pathlength be- 
tween two neighboring points on a ray to be able to 
calculate the change of optical depth along the ray 
for all wavelengths. In addition to that you must de- 
termine the angles of intersection with a given layer 
for all characteristics in order to perform angular in- 
tegrations. 

A photon in a gravitational field not only expe- 
riences a wavelength shift but also is deflected since 
it moves on a null-geodesic for the given spacetime. 
Since we are employing a characteristics based solu- 
tion we can fully account for this effect. The spatial 
derivatives in equation (O describe the geometry of 
our characteristics. They are still geodesies and can 
be described by an affine parameter. 

In order to obtain the rays, the j- part of equa- 
tion (0 has to be integrated. For the Schwarzschild 
metric it is possible to analytically describe the pho- 
ton orbits in terms of constants of the motion. In ad- 
dition, we have to relate the affine parameter to the 
distance along the characteristics. In the one dimen- 
sional case, the integration is simple and fast. Since 



the spacetime just outside the atmosphere will, in 
general, not be flat, we extend our calculation of the 
characteristics into a regime of spacetime that can be 
considered flat - in our test calculations discussed 
below, we chose a boundary of ten Schwarzschild- 
radii - to make sure that we calculate the spectrum 
from a correct set of angles that represent the imag- 
ing of the source in curved spacetime. 

We assume vacuum conditions outside the atmo- 
sphere and, therefore, the intensities will not change 
along the characteristic outside the atmosphere, of 
course except for being redshifted due to the gravi- 
tational field what is trivial to account for. 



3. The Testing Setup 

We solve the test radiative transfer problems in a 
spherical model configuration with 50 radial points 
(layers). We assume an exponential density structure 

g(r) = go exp — (7) 

'"scale 

within the atmosphere and that the gas consists only 
of a simple two-level-atom with a wavelength inde- 
pendent background continuum. 

For a given optical depth (t) grid we integrate the 
radial grid via 



where Xk - Xo&( r ) with 6 given by (0. It should 
be noted that Xk represents only the continuum ex- 
tinction coefficient. The resulting structure is not in- 
tended to be an accurate model of a neutron star at- 
mosphere. However, it has the correct spatial dimen- 
sions and we thus use it to make predictions how GR 
will affect radiative transfer in a realistic neutron star 
atmosphere. 

In Fig. [T] we plot the radial structure of the test 
atmosphere versus the optical depth grid that we 
us ed for the calculations in s ection [4] The method 
of IBaron & Hauschildtl (|2004) makes it necessary to 
provide a spatial boundary condition for the charac- 
teristics which is done via the diffusion approxima- 
tion. We generate a simple grey temperature struc - 
ture with the Hopf-function dChandrasekhaii Il950h . 
To describe the single spectral line of the two-level- 
atom, we define a wavelength as the center of 
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Fig. 1. Radius is plotted over optical depth. The op- 
tical depth was calculated from the wavelength inde- 
pendent continuous opacity. The atmosphere is about 
60 meters thick, but the layers with an optical depth 
around one lie just centimeters below the outermost 
layer. 



the line and use a Gaussian profile centered on this 
wavelength: 

y, , Wiine , /I - Aline x 2 /n\ 

OCT) = — — exp ( Y, (9) 

Wljne 

with (L)]i ne being the width of the Gaussian. We de- 
scribe the opacity associated with the line via: 

<&{X) 

^line(T,T) =^(T)/?i ine - , (10) 

J 0(T)dT 

whereby the /?ii ne factor determines the strength of 
the line relative to the continuum. 

It should be noted that the Gaussian width of 
the line is only 0.01 A. This is very small and does 
not represent a line width one would expect in neu- 
tron star atmosphere due to the large temperatures 
and pressures present in such an atmosphere . The 
small width was chosen to highlight the effects of 
GR on radiative transfer. Since the atmosphere of a 
neutron star is only a few centimeters thick, the in- 
trinsic wavelength shifts within the atmosphere are 
very small and the GR effects can be tested best with 
very narrow lines. 

A detailed treatment of radiative transfer es- 
pecially in the general relativistic environment of 
a neutron star atmosphere is desirable, since con- 
straints on the mass-radius relation are needed for 
the understanding of neutron star interiors and its 



equation of state. The constraint should be as strict as 
possible and, therefore, the radiative transfer should 
be as sophisticated as possible. 

Realistic models will need a multidimensional 
description and have characteristics in the atmo- 
sphere that extend over larger portions of space- 
time and hence have a larger intrinsic wavelength 
shift. Furthermore there may be configurations of 
blended lines that create a rapidly changing opac- 
ity as seen, e.g., in the U V spectra of classical novae 
dHauschildt etallll995l) . 

Furthermore, in realistic models of accretion 
columns on neutron stars the atmosphere extends 
over larger regions of spacetime and the intrinsic line 
shifts will be much larger. Therefore, there are phys- 
ical systems where you expect, that detailed general 
relativistic radiative transfer is important. Besides 
this method is not limited to static atmospheres but 
can also be applied to gamma-ray-bursts, accretion 
scenarios, or neutrino transport in early phases of 
core collapse supernovae. 

Since we treat the radiative transfer problem with 
scattering, we need to specify the quantities e K and 
eii ne to define the scattering albedo. The true absorp- 
tion and the scattering part of the continuum opacity 
can now be expressed through: 

k k (t) = e K XtAj) (11) 
<7>(t) = (1 - eMr) , (12) 

whereas for the line opacity we have: 
*iine(T, X) = e]Wkriine(T, X) (13) 

CHineCr, X) = (1 - eiineklineOl", /i) . (14) 

The total opacity can now be given as: 
ATtotaiCr, X) = k k {t) + ct k {t) 

+*5Bne(T, X) + cr line (T, X) (15) 

while the emissivity is: 

»7total(T, X) = (k k {t) + Kiine(T, Xj) B(T(t)) 

+(o- k (t) + 0-l ine (T, X)) J(t, X). 

(16) 

Note that the continuous opacity is constant over the 
wavelength range of interest and that the scattering 
was assumed to be coherent for simplicity. 
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Fig. 2. Results for non-scattering model line and 
continuum. 
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Fig. 3. The model line and the continuum are scatter- 
ing with ejjne = 0.01 and e K = 0.01 respectively. 



4. Results 

We calculate the emerging spectra of the model line 
for various combinations of scattering parameters for 
the continuum and for the line - hence providing 
a NLTE treatment of the line transfer. In Figs. [2] to 
[6] we always show two cases with different gravita- 
tional masses - one with M = and the other with 
M = Mq. The wavelength scale at the bottom cor- 
responds to the massless case and the upper scale 
to the one solar mass case. The emerging line pro- 
files are plotted over each other in order to be eas- 
ily compared. All calculated spectra include the full 
treatment of scattering unless it is indicated that a 
scattering parameter was set to one. 

We include the massless case in order to verify 
the code by comparing the results to the thoroughly 
tested special relativistic code and obtained the same 
results with both methods. In addition, all other tests 
such as omitting the line and recovering a flat con- 
tinuum for constant thermal sources, sudden changes 
in the wavelength resolution, and wavelength depen- 
dent boundary conditions produced correct results. 

In Fig. [2] the continuum as well as the line are 
purely thermal. In the massless case this results in a 
symmetric absorption line while in the general rela- 
tivistic case for one solar mass the line profile de- 
forms and becomes asymmetric with an extended 
wing to the red, although the effect is very small. 
The line profiles are normalized to the continuum 



and the equivalent widths of the two lines are dif- 
ferent. Although the radial structure and the run of 
the opacities is exactly the same in both cases, the 
effective opacity is different, due to the factor Aax in 
equation|2] 

In Fig. [3] the line and the continuum are scatter- 
ing with eii ne = 0.01 and e K = 0.01. In the mass- 
less case this results in a symmetric absorption pro- 
file with line wings slightly in emission, while in the 
general relativistic case the line profile is very asym- 
metric with an emission feature on the blue side and 
an extended red absorption wing. 

The model parameters used to generate Fig. |4] 
resemble those of Fig. [3] The only difference is a 
stronger scattering in the continuum e K — 1.0 x 10~ 6 . 
The strength and shape of the blue emission depends 
on the scattering in the continuum. For a strong scat- 
tering continuum even in the massless case an emis- 
sion feature appears in both wings of the li ne. It can 
be exp lained with the Schuster-mechanism (Mihalas, 
1 19781) . In the relativistic case this mechanism ap- 
pears to be amplified and strongly distorted in the 
red wing. To rule out any effect of line scattering, 
we have calculated a model line with pure scatter- 
ing (e K = 0) in the continuum and no scattering in 
the line. The result is shown in Fig. [5] The massless 
case shows the expected absorption profile with the 
wings in emission. In the relativistic case, the shape 
of the blue emission feature does not change com- 
pared to the cases with less scattering. However, its 
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Fig. 4. The model line scatters with ei; ne = 0.01 while Fig. 5. The model line is nonscattering with a com- 
the continuum scattering factor is e K — 1 .0 x 10~ 6 . pletely scattering - e K = - continuum. 



slope extends way farther into the blue, while on the 
red side of the line a very broad emission feature ex- 
tends much farther to the red, although the peak is 
not as high as the blue peak. In addition, we have 
calculated a purely scattering (eij ne = 0) model line 
combined with a purely absorptive continuum. This 
decouples the model line from the thermal pool. The 
results are shown in Fig. [6] Due to the scattering, the 
line is very strong, but the deformation of the lines 
wings is essentially the same as in Fig. [2] where the 
line did not scatter. 

The effects of GR for the given atmosphere are 
most obvious for a scattering continuum. The result- 
ing lineshape is strongly asymmetric. Scattering in 
the line has no discernible effect on the shape of the 
line. The actual deformation due to GR effects is also 
very small for pure linescattering. 

The emergent lineshapes are related in princi- 
ple to P-Cygni profiles, since there is also a con- 
stant wavelength shift - the dopplershift - through- 
out the atmosphere. However there is a fundamental 
difference between the two cases, since in an expand- 
ing atmosphere where P-Cygni profiles are observed 
there is a point of last contact of a photon with the at- 
mosphere, hence the information of the dopplershift 
at the point of emission in respect to the observer is 
conserved. In the gravitational field the wavelength 
of the photon gets shifted even without any interac- 
tion with the atmosphere. 
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Fig. 6. The model line is completely scattering - 
eiine = - while the continuum is not scattering at 
all. 



5. Grey continuum models 

As another application to general relativistic transfer 
we calculated grey continuum models. To do so we 
changed our wavelength resolution and omitted the 
model line from our calculations CRii ne = 0). In the 
following we present the emerging continuous spec- 
tra for varying values of e K . 

The temperature structure of the model is grey 
with an effective temperature of 10 4 K (the absolute 
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value of the effective temperature is not important for 
the testing). The emergent spectra are nearly black- 
body spectra with only a slightly distorted overall 
shape as they are a little bit broader than the black- 
body. 

To determine the effective temperature from an 
observation one would have to correct for the grav- 
itational redshift and then try to fit a blackbody to 
the spectrum. This procedure will only deliver phys- 
ically relevant results if the scattering conditions in 
the atmosphere are known. As long as all photons 
are absorbed and none scattered - e K — 1 - (see 
Fig. Q the emergent spectrum (corrected for grav- 
itational redshift) is fitted with a blackbody with a 
temperature that is the same as the effective temper- 
ature of the model. For models with non-zero scat- 
tering - e K = 0.1 in Fig. [8] and e K - 0.001 in Fig. 
[TOl - the emerging spectra are much bluer than the 
blackbody fit with the model effective temperature. 
In Fig. [9] and Fig.QT]we corrected the temperatures 
of the blackbody fits to match the emerging spectra. 
The apparent temperatures are much higher than the 
effective temperature of the model. 

To prove the consistency of our results, we plot 
the thermalisation depth t± - the optical depth where 
the J,i = - over the scattering parameter e K . Since 
the emergent spectrum is Planckian, r t h is the opti- 
cal depth where the temperature of the atmosphere 
equals the temperature of the blackbody fit. The re- 
sults are plotted in Fig. [12] 

The results are consistent with the simple model 
r t h = A=. Hence, in the GR case the determination 
of effective temperatures via blackbody fits or sim- 
ple radiative transfer models that neglect scattering 
will be inconsistent and result in systematic errors 
of the effective temperatures. Therefore, it is desir- 
able to solve the scattering problem in a GR envi- 
ronment selfconsistently in a full blown atmospheric 
code with the method of solution that we present in 
this paper. 

6. Conclusion 

We have developed a method to solve continuum and 
line radiative transfer problems in spherically sym- 
metric spacetimes that fully accounts for general rel- 
ativistic effects and can account for scattering in the 
continuum and the lines. It uses a comoving frame 
wavelength formalism that allows resolution of spec- 
tral lines throughout the atmosphere without signif- 



Blockbody: 10000 K 




X [A] 

Fig. 7. The spectrum was corrected for the gravita- 
tional redshift and fitted with a blackbody with the 
known effective temperature of the model - T e g = 
10 4 K. The continuum has e K — 1 The blackbody fits 
the model reasonably well. 



icantly increasing the number of wavelength points. 
The method was developed and tested in static neu- 
tron star-like atmospheres, but is generally applica- 
ble to general relativistic systems. Only the photon 
orbits and the wavelength coupling term ax would 
be different in other GR systems. The test models 
provide an illustration of possible results for realis- 
tic model atmospheres for neutron stars. The results 
show that the emergent line profiles in general rel- 
ativistic atmospheres cannot be described in detail 
with non-relativistic radiative transfer. The influence 
of the continuous scattering opacity on the shape of 
the lines is large. 

The apparent effective temperature of continuous 
spectra also depends strongly on the strength of the 
scattering. Therefore, it is necessary to include the 
treatment of scattering in the radiative transfer solu- 
tion in order to obtain a consistent physical model 
of a neutron star atmosphere and similar cases. The 
method that we have presented here is a first step 
to develop a thorough treatment of general relativis- 
tic atmosphere models in 3D. It can be directly ap- 
plied to multi-level NLTE calculations of relativistic 
neutron star atmospheres, which we will present in a 
subsequent paper. 
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Fig. 8. The spectrum was corrected for the gravita- 
tional redshift and fitted with a blackbody with the 
known effective temperature of the model - T e ff = 
10 4 K. The continuum has e K = 0.1 The blackbody 
doesn't fit the model spectrum. The apparent tem- 
perature is higher than the model temperature. 



Fig. 10. The spectrum was corrected for the gravita- 
tional redshift and fitted with a blackbody with the 
known effective temperature of the model - T e ff = 
10 4 K. The continuum has e K — 0.001 The blackbody 
doesn't fit the model spectrum. The apparent temper- 
ature is higher than the model temperature. 
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Fig. 9. The model spectrum is the same as in Figure 
[8] but this time the temperature of the blackbody was 
chosen to fit the spectrum. 



Blackbody: 21400 K 
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Fig. 11. The model spectrum is the same as in Figure 
[TOl but this time the temperature of the blackbody 
was chosen to fit the spectrum. 
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